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Preface 
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Abstract 

i 

The  report  dealt  with  a  tradeoff  between  time- to- turn 
and  the  change  in  specific  energy  during  a  turn  of  180°, 
for  an  aircraft  of  high  specific  power.  This  type  of 
aircraft  possesses  the  capability  to  sustain  flight  at 
the  corner  velocity  where  the  aircraft  flies  at  maximum 
lift  coefficient  and  maximum  load  factor  simultaneously. 
However ,  the  classical  necessary  conditions  breakdown  on 
this  corner  velocity  arc  and  an  addition  constraint  set 
must  be  defined  to  determine  the  optimal  control 
histories.  The  report  first  defines  the  necessary 
conditions  for  a  generic  optimal  control  formulation  with 
two  state-dependent  inequality  constraints  and  then 
applies  the  formulation  to  a  high  specific  power  aircraft 
problem.  The  result  is  a  three-point  boundary  value 
problem  including  a  discontinuous  interior  corner  time  at 
the  beginning  of  sustained  corner  velocity  flight.  All 
state  and  costate  derivatives  and  end  conditions  are 
presented  with  numerical  methods  for  determining  minimum 
time#  maximum  energy  and  maximum  energy  gradient 
solutions.  Anticipated  results  are  provided  for  the 
formulations  beginning  at  initial  conditions  both  above 
and  below  the  corner  velocity.  Results  are  presented  for 
the  minimum  time  solution  only.  Recommendations  for 
further  work  in  this  same  area  are  also  provided. 
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HIGH  SPECIFIC  POWER  AIRCRAFT  TURN 
MANEUVERS:  TRADEOFF  OF  TIME -TO- TURN 
VERSUS  CHANGE  IN  SPECIFIC  ENERGY 

I.  INTRODUCTION 

Since  World  War  I*  great  strides  have  been  made  in 
the  field  of  energy  management  of  combat  aircraft.  Most 
of  the  presently  used  tactics  were  developed  during  the 
early  portion  of  this  century*  but  their  detailed 
explanation  and  simulation  has  only  been  accomplished  in 
the  past  25  years.  As  more  progress  is  made  in  aircraft 
design  and  capabilities*  new  avenues  of  research  in 
energy  management  are  opened  for  examination. 

Much  work  has  been  done  in  the  areas  of  fuel  and 
energy  management  for  fighter  aircraft  of  the  1960's  thru 
1970's  vintage*  particularly  the  area  of  minimum  time-to- 
turn.  In  research  of  any  particular  aircraft's  minimum¬ 
time  turn  performance*  a  solution  is  sought  whereby  the 
angle  of  attack*  bank  angle  and  throttle  histories  create 
an  exchange  of  potential  and  kinetic  energy  to  produce 
the  highest  turn  rate  possible  at  all  points  in  time. 
For  every  aircraft*  there  is  a  velocity  which  is  a 
function  of  altitude*  at  which  the  maximum  possible  turn 
rate  occurs.  It  is  that  airspeed  at  which  maximum  lift 
coefficient  and  maximum  load  factor  occur  simultaneously 
and  it  is  known  as  the  corner  velocity.  Beginning  at  a 
particular  initial  airspeed  and  altitude*  flying  as  close 
as  possible  to  this  velocity  at  all  times  gives  the 
minimum  time  turn  for  that  aircraft. 
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Historically#  aircraft  have  possessed  excessive 
weight  and  induced  drag  and  insufficient  thrust  to 
sustain  flight  at  their  corner  velocities  although  they 
can  instantaneously  pass  through  this  condition  during  a 
maximum  performance  maneuver.  However#  as  technology 
advances#  lighter#  more  aerodynamically  efficient 
aircraft  with  increased  thrust  are  developed  which  can 
more  closely  approximate  sustained  corner  velocity 
flight. 

In  the  near  future#  fighter  aircraft  may  produce  more 
than  the  minimum  required  specific  power  to  sustain 
corner  velocity  flight  and  will  be  able  to  achieve  the 
maximum  possible  turning  rate  at  less  than  maximum 
throttle  setting.  Should  this  capability  exist  for  a 
particular  aircraft#  it  would  be  of  interest  to  examine 
the  tradeoff  between  time  allotted  to  complete  a  180° 
velocity  vector  heading  change  and  the  resultant  specific 
energy  change  at  the  end  of  the  turn#  for  turns  that 
begin  both  below  and  above  the  corner  velocity.  The 
examination  of  this  tradeoff  and  the  set  of  optimal 
control  histories  which  produce  it  is  the  objective  of 
this  thesis. 

This  work  is  an  extension  of  the  work  done  by 
Humpureys  (Ref  5)  and  Anderson  (Ref  1).  In  this  effort# 
the  aircraft  model  of  Humphreys  and  the  kinematic  optimal 
control  equations  developed  by  Anderson  are  blended  with 
desired  initial  and  final  flight  conditions  to  form  a 
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pair  of  Threa  Point  Boundary  Value  Problems  (TPBVP);  one 
for  minimum  time  (free  final  time)  and  one  for  maximum 
energy  gain  (fixed  final  time).  Prom  these  formulations, 
the  tradeoff  is  determined. 
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II.  Statement  of  the  Problem 
Connection  of  Energy  and  Time  in  a  Turn 

As  stated  in  the  introduction,  the  objective  of  this 
thesis  is  to  examine  the  tradeoff  between  time  to 
complete  a  180°  velocity  vector  heading  change  and  the 
resultant  specific  energy  change  from  the  beginning  of 
the  turn  to  the  end.  Therefore,  the  tradeoff  criterion 
between  energy  and  time  must  be  defined. 

The  total  energy  of  an  aircraft  is  the  sum  of  its 
potential  and  kinetic  energy  components  and  can  be  stated 
as: 


E*r  »  mgh  +  %mV2  (2-1) 

which  has  units  of  foot-pounds.  The  specific  energy,  or 
the  energy  altitude  of  an  aircraft  is  the  total  energy 
per  pound  of  aircraft  weight  and  can  be  expressed  as: 

E  -  h  +  Y2  (2-2) 

2g 

which  has  units  of  feet.  If  an  aircraft  begins  a  turn 
with  some  initial  specific  energy,  E},  the  effects  of 
climbing  and  diving,  drag  and  thrust  application  will 
transition  the  initial  aircraft  energy  state  to  a  final 
energy  state,  Ef,  in  a  time,  tf. 

For  a  minimum  time  turn,  energy  is  sacrificed  for  the 
sake  of  continuous  maximum  turn  rate.  In  this  case,  the 
final  energy  Emin  will  be  the  lowest  final  energy  of  any 
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of  the  turns  to  bs  examined  and  it  will  occur  at  a  time 
tnin*  tha  minimum  tima  to  turn  to  tha  final  conditions. 
Zf  tha  tima  allotted  to  turn  is  increased  to  tf  >  tmin, 
than  tha  opportunity  exists  to  perform  tha  same  turn  in  a 
more  energy  conscious  manner  and  therefore#  Ef  >  Emin. 

The  more  time  allotted  for  the  turn#  the  more  energy  can 
be  conserved  or  gained. 

Definition  of  the  Tradeoff 

Xt  is  postulated  that  a  final  specific  energy#  Etrade' 
occurring  at  a  final  time  tirade  will  maximize  the 
specific  energy  gradient  defined  as: 

EG  »  E* "  Ei  (2-3) 

Etrade'  ttrade  and  ZGmax  be  graphically  or 

numerically  determined  by  a  line#  eminating  from  E*  at 
t*0  being  placed  tangent  to  a  final  specific  energy 
versus  final  time  plot.  This  is  qualitatively  shown  in 
Fig  2-1.  An  additional  method  for  determining  EGmax  is 
to  plot  the  energy  gradient  versus  final  time  and  pick 
the  maximum  value.  This  can  be  seen  from  Fig  2-2. 

The  determination  of  ttrade'  Etrade  *nd  therefore 
E6max  for  turns  beginning  at  equal  specific  energy 
states#  both  below  and  above  the  corner  velocity#  is  the 
goal  of  this  effort.  This  could  be  done  using  a  single 
problem  formulation#  such  as  maximize  EG#  instead  of 


Final  Specific  Energy,  E  (feet) 


•ade 


Figure  2- 


Time,  t  (sec) 


.  Specific  Energy  Gradient  vs.  Time-to-Turn 
depicting  EGmax. 
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>$•/,  defining  a  minimum  time  turn  and  several  maximum  energy 

turns,  but  it  would  not  be  as  useful  a  result  for 
comparative  purposes  as  the  formulation  chosen.  This 
method  unambiguously  determines  the  time-energy  tradeoff. 
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III.  Generic  Optimal  Control  Formulation 
The  system  to  be  optimally  controlled  will  be 
represented  by  the  dynamic  differential  system  (Ref  1) 

x  ■  f(x,u,t)  (3-1) 

with  initial  and  final  time  conditions 

x(t0)  -  x0  ,  K [x( tf ) #  tf ]  -  0  (3-2) 

Both  x  and  x  are  n-dimensional  state  vectors,  u  is  an  m- 
dimensional  control  vector,  K  is  an  s-dimensional  vector 
of  final  conditions  and  t  represents  time.  The  problem 
objective  is  to  transfer  the  system  from  the  initial 
conditions  to  the  final  conditions  while  minimizing  the 
Bolza  form  payoff  function 

J  ■  ♦  V1k[,(t(),t(] 

+  fj*  [F(x,u,t)  +  AT(f-x)]dt  (3-3) 

G  is  a  scalar  function  of  final  time  and  final  time  state 
conditions  and  F  is  a  scalar  integral  function  subject  to 
state  and  control  values  and  time.  The  vector  uT  is  an 
s-dimensional  constant  Lagrange  multiplier  vector,  and 
\T  is  an  n-dimenional  Lagrange  multiplier  costate  vector 
otherwise  known  as  the  influence  function  vector.  The 
Hamiltonian  is  expressed  as: 

H(  x,  u,  t)  »  F(x,u,  t)  +  ATf(x,u,t)  (3-4) 
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The  control  vector,  u,  is  subject  to  both  non-state- 
dependent  and  state-dependent  inequality  constraints. 
The  non-state  dependent  constraints  do  not  affect  the 
control  formulation,  only  the  control  vector  bounds. 
They  will  be  considered  separately  in  section  4.  For  the 
class  of  problems  discussed  here,  it  is  postulated  that 
•the  ifch  element  of  the  control  vector  is  subject  to  two 
state-dependent  inequality  constraints 

Cl(x,Ui,t)  £  0,  C2(x,ui,t)  4  0  (3-5) 

Adjoining  these  constraints  to  the  Hamiltonian,  it 
becomes 

H  =  F(  x,  u,  t )  +  ATf(x,u,t)  +  iC]_  ( x,  u£ ,  t) 

+  02C2(x,ui,t)  (3-6) 

where  m  and  yu2  are  scalar  Lagrange  multipliers.  The 
conditions  on  m  and  fJ2  are 

/lj  ■  0  if  Cj(x,ui,t)  £  0 

.  1-1,2  (3-7) 

/ij  ^  0  if  Cj  ( x,  ui ,  t)  «  0 

if,  and  only  if,  both  and  C2  are  not  both  zero 

simultaneously.  If  only  one  constraint  is  zero  at  any 
time,  applying  the  Euler-Lagrange  optimality  condition, 

Hu  ■  0  (3-8) 

produces  the  impulse  response  functions  and  the  non-zero 
value  of  Jlj  can  be  found  from  the  relationship 
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v*A  Hu£  ■  Fuj  +  XTfUi  +  /ijCjm  ■  0  j  ■  1,2  (3-9) 

However,  if  both  state-dependent  inequality  constraints 
are  zero  simultaneously,  Eq  (3-9)  becomes 

Hu*  ■  Fui  +  XTfUj[  +  Z^lCiUi  +  A*2C2U1  *  0  (3-10) 

as  Eq  (3-10)  is  the  only  equation  which  yields  values  for 
Ml  or  M2*  unique  solutions  for  Ml  and  M2  cannot  be 
obtained  with  one  equation  and  two  unknowns,  and 
alternate  solution  must  be  sought. 

The  solution  selected  is  to  form  a  new  state 
dependent  equality  constraint  from  both  Ci  and  C2 
simultaneously  equaling  zero  (Ref  1:182,  2:2247,  3:117). 
This  new  constraint  has  the  form 

S(X't)  -  0  (3-11) 

However,  this  constraint  is  not  explicitly  a  function  of 
the  control  vector  and  the  optimal  controls  cannot  be 
found  from  this  expression.  Using  the  techniques  of  (Ref 
3:118),  it  is  found  that  the  qth  time  derivative  of  Eq 
(3-11)  produces  a  constraint  containing  elements  of  the 
control  vector.  This  new  constraint 

s(3)(x,u,t)  ■  0  (3-12) 

is  known  as  a  singular  arc  and  may  now  be  adjoined  to  Eq 
(3-6)  to  form  the  new  Hamiltonian 
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H  ■  F(x,u,  t)  +  ATf(x,u,t)  + /ilCl(xf  Ui,  t) 

+/U2C2(x,ui,  t)  +  fJI  3S<<3)  (x,u,  t)  (3-13) 

For  this  new  formulation,  it  should  be  understood  that 
Eqs  (3-7)  are  still  valid  if  only  one  constraint  is  zero 
at  a  time.  If  both  Ci  and  C2  are  zero  simultaneously, 
then  Eq  (3-11)  is  satisfied  and  Hi  and  n 2  are  defined  to 
both  be  zero.  The  unique  value  for  /J 3  is  determined  by 
applying  Eq  (3-8)  to  Eq  (3-13)  which  produces 

Hu  «  Fu  +  ATfu  +  Ji3Su(<3)  «  0  (3-14) 

In  addition  to  the  singular  arc  constraint,  q 

interior  point  constraints  have  been  generated  which  must 
be  adhered  to  at  the  beginning  and  end  of  the  arc,  and 

are  represented  in  vector  form  as 

M  *  [  S,  S,  ....  ]  -  0  (3-15) 

The  Hamiltonian,  H,  and  the  costate  variables,  AT  are 
discontinuous  at  the  beginning  of  the  arc,  but  continuous 
at  the  end  (Ref  2:2248).  If  the  discontinuity,  or 

corner,  happens  at  time  tc,  then  the  shifts  in  the  values 
of  H  and  AT  are  the  differences  between  their  values  just 
before  and  after  tc  and  can  be  expressed  as 

H(tc-)  „  H ( tc+ )  -  77TMt  t  *  fcc  (3-16) 

AT(tc-)  =  AT(tc+)  +  77tMx  t  -  tc  (3-17) 

where  *7T  is  a  q-dimensional  constant  Lagrange  multiplier 
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vector.  The  values  of  are  found  by  simultaneously 
solving  Eg  (3-16)  and  Eq  (3-17). 

In  summary,  given  the  system  and  conditions  of  Eqs 
(3-1)  and  (3-2),  the  payoff  function  of  Eq  (3-3)  becomes 

J  -  G  +  utK 

[P+  +  /*T(ci,C2,S<q>)]  dt  (3-17) 

and  the  Hamiltonian  is  expressed  by  Eq  (3-12).  The 

optimal  controls  are  found  by  applying  Eg  (3-8)  to  Eq  (3- 

13)  subject  to  the  constraints  of  Eqs  (3-7),  those 

following  Eq  (3-13),  and  insuring  satisfaction  of  the 

> 

Legendre-Clebsch  condition 

Hju>0  (3-18) 

The  costate  differential  equations  equal  -Hx  and  are 

At  -  -<F*  +  ATf*  +/»lClx  +M2C2x  +<l3Sx<1>)  (3-19) 

Application  of  the  transversality  conditions  yield  the 
final  time  values  of  the  costates  and  Hamiltonian  to  be 

*T(tf)  «  Gx[x(tf),tf]+  oTKx[x(tf),tf  (3-20) 

H(  tf )  -  -Gt[x(tf),tf]-  VTKt[x{tf),tf  ]  (3-21) 

which  completes  formation  of  the  three  point  boundary 
value  problem  (TPBVP).  Eq  (3-21)  applies  only  to  the 
free  final  time  problem.  The  interior  point  constraints 
of  Eq  (3-15)  and  the  corner  discontinuities  of  Eq  (3-16) 
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^  an<^  (3-17)  must  also  be  satisfied  when  switching  onto 

or  off  of  the  singular  arc  of  Eq  (3-12)  at  the  corner 
time,  tG. 
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IV.  Specific  Optimal  Control  Formulation 


Point  Mass  Aircraft  System  Definition 

To  satisfy  Eq  (3-1 ),  the  aircraft  for  both  minimum 
time  and  maximum  final  energy  turn  problems  will  be 
represented  by  the  following  point  mass  equations  of 
motion  (Ref  12) 


X  ■  V  cosy  cos 

Y  ■  v  cosy  sin  0 
h  ■  V  siny 

V  »  g|l7rzD  -  sin yj 

y  *  3L  [li  cos  0  -  cosy  ] 

v[w  J 

ip  *  -3,.  [  L  sin0  ] 
vcosy[  w  j 


(4-la) 
(4-lb) 
(4-lc) 
( 4-ld ) 


(4-le) 


(4-lf ) 


The  variables  x,  y  and  h  form  the  three  axis  position,  V 
is  velocity,  y  is  the  flight  path  angle  and  0  is  the 
heading  angle.  T  is  the  constant  maximum  available 
thrust,  D  is  drag,  L  is  lift,  W  is  a  constant  weight  and 
g  is  the  acceleration  of  gravity.  The  expressions  for 
lift  and  drag  are 


L  »  h  P0  a  V2ACL 
D  -  H  P0<rV2A(CD  +kCL2) 


(4-2a) 

(4-2b) 


where  A  is  the  reference  wing  area,  CL  is  the  lift 
coefficient,  Cdq  is  the  zero  lift  drag  coefficient  and  k 
is  the  induced  drag  parameter.  The  parameter  f)0a  is 
the  exponential  air  density  expression 
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PqO  -  p0e  8  (h-ho) 


(4-3a) 


where  $  is  a  constant  and  h^  is  the  minimum  altitude  of 
interest.  Least  squares  curve  fits  to  data  from 
atmospheric  tables  (Ref  8)  produce 


Po 

Po  o  3 


0. 0023769e 
0. 00072736e 


*  m 


30477 
h-36089 


20953 


(4“3b) 
( 4 — 3  c ) 


Where  Eqs  (4-3b)  and  (4-3c)  represent  the  troposphere  and 
stratosphere  respectively.  For  this  exercise,  only  Eq 
( 4-3 b)  will  be  employed  due  to  the  altitudes  at  which  the 
aircraft  maneuvers  will  be  examined. 

The  control  vector  is 


u  »  [  0  CL  IT  ]  (4-4) 

where  0  is  the  bank  angle,  is  the  lift  coefficient 
and  7 T  is  the  throttle  setting.  The  constraints  on  the 
control  vector  are 


-180°£  $  5+  180° 

(4-5a) 

Cl  -  (CL  -  CLmax)  £  0 

(4-5b) 

C2  *  (C&  -  2WNmax)  0 

(4-5c) 

P0V2A  a 

c3  -  CL>0 

( 4-5d ) 

-0.5  £  ir  z  l.o 

(4-5e) 

where 
and  N 


Ckmax  *s  max^jnum  aerodynamic  lift  coefficient 
max  is  the  maximum  allowable  normal  load  factor. 
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From  Eq  (4-3)  it  is  seen  that  the  bank  angle  is 
unbounded,  the  lift  coefficient  muut  always  be  greater 
than  zero,  less  than  both  the  maximum  aerodynamic  lift 
coefficient  and  that  for  achieving  maximum  load  factor 
and  the  throttle  setting  may  vary  from  half  reverse 
thrust  to  full  forward  thrust  settings. 

High  Specific  Power  Aircraft  -  Corner  Velocity  Constraint 
If  both  Eqs  ( 4-5b, c)  are  simultaneously  zero,  Eq  (3- 
11)  is  formed  and  is 

S(x,t)  3  aV2  -  2WNmax  *  0  (4-6) 

PoAcLmax 

This  expression  must  be  satisfied  if  sustained  maximum 
turn  rate  flight  at  the  corner  velocity  is  to  be  part  of 
the  optimal  trajectory.  However,  it  is  readily  seen  that 
Eq  (4-6)  is  not  explicitly  a  function  of  any  control  and 
therefore  successive  time  derivatives  must  be  taken  until 
Eq  (3-12)  is  formed.  This  equation  is 

S(x,  t)  ■  V2  Chh  +  2  aW  ■  0  (4-7) 

Substitution  of  Eq  (4-lc,d)  into  Eq  (4-7)  and  rearranging 
the  resultant  expression  yields 

su.u.ti  .  •[!$*-  «*«r 

This  expression  is  a  function  of  throttle  setting  and 
lift  coefficient  and  satisfies  the  requirements  of  Eq  (3- 


(1 


-tv2) 

75 


m 


(4-8) 
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12).  Therefore,  only  one  interior  point  constraint 
exists,  Eq  (4-6),  and  Eq  (3-15)  becomes 

M  »  S(x,t)  -  <rv2  -  2WNmax  -  0  (4-9) 

Po*cLroax 

The  Hamiltonian  of  Eq  (3-13)  can  now  be  formed  and  is 
H  *  AgX  +  XyY  +  Xfch  +  X  y  V  4*  X y  y  ^  X^  ^ 

+  -  IW>+  «39[2^£  -  .iny(l-^v2>]  <4*10> 

2g 

For  this  type  of  aircraft  problem,  there  is  no  integral 
function  F(x,u,t)  and  the  constraint  Ci  of  Eq  (4-5b)  is 
not  state  dependent  and  need  not  be  adjoined  to  the 
Hamiltonian.  Neither  M  nor  H  is  explicitly  a  function  of 
time  and  therefore  from  Eq  (3-16)  it  is  seen  that 

Mt  |  t-tc  *  0  H(tc”)  -  H ( tc+  j  (4-11) 

Ht  ■>  0  .*.  H(tc)  ■  H(tc)  -  H(tf)  *  constant  (4-12) 

Recalling  Eq  (3-17)  and  applying  Eq  (4-9),  it  is  seen 
that  only  Xh  and  Xv  will  have  interior  point 
discontinuities  which  are 

*h<tc  )  *  *h<tc+)  +  rifiav2  (4-13a) 

^-v^c  )  ■  Xv(  tc+)  +  2W<rV  (4-13b) 

Upon  substitution  of  Eq  (4-13)  into  Eq  (4-11)  it  is  found 
that 
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(4-14) 


Resubstitution  of  Eq  (4-14)  into  Eq  (4-13)  produces 

Ah(  fcc")  *  Ah(  tc+)  ”  Av(  tc+)®^  l 

2 


(4-15a) 


Xv(te)  ■  0 


(4-15b) 


The  discontinuity  in  Av  and  the  dependence  of  Ah(tc“)  on 
Av(tc*)  makes  forward  integration  of  the  differential 
equations  impossible  and  backward  integration  from  the 
final  conditions  to  the  initial  conditions  is  required. 
This  will  be  discussed  in  Sec  5. 

Optimal  Controls 

The  optimal  controls  are  defined  by  applying  the 
Eulor-Lagrange  condition,  Eq.(3-8)  and  the  Legendre  - 
Clebsch  condition,  Eq  (3-18)  to  the  Hamiltonian  of  Eq  (4- 
10).  To  find  the  optimal  bank  angle,  0  ,  two  equations 
must  be  satisfied.  They  are 

H*  =  A y3.  L(-Sin  0  )  +  g  \  Lcos 0  ■  0  (4-16a) 

9  Yv  w  vcosy  ♦w 


H4  ■  -  Ay  sin  0  +  A#£OS^  =  0 

cosy 

H*  -  -  Ay  cos 0  -  A*  sinft  >  0 

cosy 


(4-16b) 


(4-17) 


In  order  for  both  of  these  conditions  to  be  satisfied 
simultaneously,  the  only  solution  is 


cos  0 


,2  JJuLf 

Kr  lcos  / 


,  sin0  ■ 


-  \*/cos 


2  (4-18) 


19 


AP I T/GAE/AA/7  8  D- 5 


Upon  substitution  of  the  non-uniqut  cos0  and  sin0 


expressions  into  Eq  (4-17)  it  is  seen  that 


\ /  Xy2  + 


cosy 


(4-19) 


thus  insuring  a  minimizing  control  at  all  times.  Eq.  (4- 


18)  will  produce  the  optimal  bank  angle  history. 


To  define  the  optimal  lift  coefficient,  again  two 


conditions  must  be  adhered  to.  They  are 


-aXvDCL  +  2_AyLcLcos0  +  Jk  LcLsin0 


HCr  =  -axvucL  + 

u  w  wv 


WVcos 


+  a  <'2LCL  -  a  l<3DcL  -  0  (4-20) 

HcLCL  *  -  a  ^vDClCl  ‘  g  ''3DCLCt  ~  0  <4‘21) 


Rearranging  and  substituting  Eq  (4-17)  yields 


«CL  -  -(Xv  +  ^3)°CL  ”  ~  U2)^CL  "  0  (4-22) 


HCLCL  *  -(  Ay  +  ^3)dClCl  £  0 


(4-23) 


For  an  interior  Cl  control  which  is  within  the 


constraints  of  Eqs  (4-5b,c,d)  both  fi 2  and  H 3  are  zero. 


Therefore,  Eq  (4-23)  states  that  for  an  optimum  interior 


control 


Av  £  0 


(4-24) 


as  DClCl  *s  alway®  positive.  Eq  (4-22)  produces 


optimum  interior  control 
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cL 


zHjsl 

2kvAv 


(4-25) 


If  any  of  the  bounds  of  Eq  (4-5b,c,d)  are  exceeded  by  Eq 
(4-25),  0^  is  set  to  the  boundary  value  exceeded.  If 
maximum  load  factor  flight  occurs,  Eq  (4-5c)  is  zero,  and 
the  non-zero  value  of  M 2  calculated  from  Eq  (4-22)  to 
be 


H  2  - 


4kWNmax^v  + 
p0  a  v2a  T 


(4-26) 


To  determine  the  condition  of  HcL  during  corner  velocity 
flight,  the  throttle  history  must  be  examined  first. 

The  two  final  equations  which  must  be  investigated 
for  the  optimal  thrust  control  are 

H „  *  (  Xv  +  V3)ql  -  0  (4-27) 

W 

and 

Hvy  =  0  (4-28) 

Eqs  (4-27)  and  (4-28)  demonstrate  that  the  throttle  is  a 
singular  control  because  it  does  not  appear  in  Eq  (4-27). 
Therefore,  unless  // 3  equals  \v,  which  can  only  happen 
during  sustained  corner  velocity  flight,  the  throttle 
setting  cannot  be  an  interior  value,  but  must  be  a 
minimum  or  maximum  boundary  value  depending  upon  the 
value  of  Ay  Using  Eq  (4-5e)  and  Eq(4-27)  it  is 


determined  that 
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*  \ 
V  r 


ir  ■  1.0  for 

Xv  «  0 

(4-29a) 

7T  •  -0.5  for 

Xv  >  0 

( 4-29b) 

If 

Xv  equals  zero, 

which  happens  as  the  corner 

velocity 

arc 

is  entered  as  shown  by  Eq  (4-15b),  then  the 

throttle 

value  is  determined  to  be 

If  ■  1.0  for 

Xv  *£  0 

(4-30a) 

7 T*  -0.5  for 

Xv  *  0 

( 4  — 3  0  b ) 

For  flight  along  the  corner  velocity  arc,  the  partial  or 
interior  throttle  setting  can  be  determined  from  Eq  (4-8) 
to  be 

V  -  i  [Dmax  +  W  siny  (1-  *V2}  1  (4-31) 

T  1  2g 

For  this  interior  control  to  be  a  viable  solution,  from 
Eq  (4-27),  it  is  seen  that 

/i3  =  -  Xv  (4-32) 

Upon  substitution  of  Eq  (4-32)  into  Eq  (4-22),  it  is 
determined  that  along  the  corner  velocity  arc 

HCl  -  "  MfsLLcL  £  0  (4-33) 


as  all  of  the  components  of  HcL  are  always  positive. 
Therefore,  Eq  (4-33)  demands  that  be  the  optimum  Cl 
for  this  condition  which  corroborates  the  requirements  of 
Eq  (4-6). 
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Costa ts  Differential  Equations 

Of  the  three  Euler-Lagrange  conditions,  two  have  been 
developed)  one  for  the  state  derivatives  of  Eqs  (3-1)  and 
(4-1)  and  a  second,  Eg  (3-8),  to  determine  the  minimising 
optimal  controls  just  discussed.  The  third  condition,  Eg 
(3*19),  deals  with  the  costate  differential  system,  and 
produces  the  following  eguationsi 


**  0 
Xy  *  0 

*h  -  |  [°h(  XV  +  A<3>+I*h(&l.  “  **2>] 

Xv  -  -3  f  Dv(  Xv  +  /i3)+Lv(SM.  -  M2)\ 

W  1  V  J 

+  I(  Xy  Y  +  \A>)-9in  y  (  Xh  +  M3  V) 
V 

-  cos  y  (  xx  cos \ft  +  Xy  sin^  ) 


siny  |  V(Xxcos^  +  Xy  sin^)  -  Xy  3. 

cosy  g  [  xv  +  M3(l-  0v2)]  -  Xhv 

3 


X^tan  y 


(4-34) 

(4-35) 

(4-36) 

(4-37) 


(4-38) 


cos  y  (  X  x  sin  iff  - 


cos  ip  ) 


(4-39) 


Aircraft  Model 

The  coefficients  used  to  complete  the  point  mass 
aircraft  model  are  identical  to  those  of  (Ref  5:97).  The 
quantities  in  table  (4-1)  were  used  to  calculate  the 
initial  altitude  and  airspeed  values  of  Table  (4-2). 
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Name 

parameter 

value 

Thrust#  max 

T 

18,225  lbs* 

Weight 

W 

12,150  lbs 

Maximum  load 
factor 

Nmax 

7.22 

Reference 

Wing  Area 

A 

237  ft2 

zero  lift 
drag  coefficient 

n 

o 

0 

0.02 

induced  drag 
coefficient 

k 

0.05 

maximum 

lift  coefficient 

Cl*max 

1.0 

*  -  This  thrust  corresponds  to  T/W  *1.5 


Table  4-1.  Point  Mass  Aircraft  Model  Parameters 
Transversality  Conditions 

As  stated  in  Sec.  3,  the  transversality  conditions# 
Eqs  (3-20)  and  (3-21 )#  complete  the  formation  of  the 
TPBVP.  For  the  minimum  time  and  maximum  energy  turns  to 
be  investigated#  the  state  dependent  conditions  of  Eq  (3- 
2)  are  as  shown  in  Table  (4-2).  From  that  table#  it  is 
evident  that  for  all  cases  examined. 


|  X( tf ) #  tf  )  J ■ 

y 

tjt-ir 

(4-40) 


t*tf 


where  W  in  this  case  represents  180°  in  radians. 
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Table  4-2.  Initial  and  Final  Conditions  for  Minimum  Time  and  Maximum  Energy  Problems 
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Employing  Eq  (3-20),  it  is  determined  that 

[  Xx  Xy  Xh  Xv  Xr  X#J  t*tf  ■  [  0  0  flh  Gv  ^  *■'2  ]  (4-41) 

which  are  not  a  function  of  the  state  final  conditions, 
K.  This  result  holds  for  all  three  problem  formulations. 
It  is  noted  that  from  Eqs  (4-34),  (4-35)  and  (4-41)  that 
Xx  and  Xy  are  always  zero,  have  no  effect  upon  the 

problem  solutions  and  may  be  eliminated  from  the 
formulations  from  here  on.  For  simplification,  Eq  (4-39) 
now  becomes 

X*  -  0  (4-42) 

Eq  (4-41)  shows  that  values  for  X^  and  Xv  at  the 

final  time  can  be  explicitly  calculated,  but  values  for 
XY  and  X*  must  be  judiciously  selected.  No 

additional  help  can  be  provided  for  selecting  initial 
guesses  for  the  finals  values  of  Xy  and  X*  for  a 

fixed  final  time  problem;  only  that  they  together  must 
generate  a  feasible  final  time  bank  angle. 

However,  for  both  free  final  time  problems,  using  Eqs 
(3-13),  (3-21),  (4-8),  (4-40)  and  (4-41),  it  can  be  shown 
that  at  t  ■  tf 

Xyy  +  X*(//  +  Gtf  -  0  (4-43) 

As  noted  in  (Ref  1*183),  when  a  singular  arc  of  this 
type  is  part  of  the  optimal  trajectory,  the  trajectory 
terminates  on  the  arc.  This  is  helpful  for  two  reasons: 
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1)  the  final  velocity  ia  the  corner  velocity  and  2)  the 
final  load  factor  is  maximum  allowable  load  factor*  Upon 
substitution  of  these  two  facts  plus  Eqs  (4-le,f),  (4-17) 
(4-19)  and  quantities  from  table  (4-2),  it  is  seen  that 

Jxr2  +  A,2  +M  -  St,  -  0  <«-♦«> 

With  one  equation  and  two  unknowns,  one  variable  must 
be  selected.  To  help  make  the  best  selection,  two 
observations  can  be  made;  1)  based  on  Eq  (4-40)  a 
positive,  right  turn  must  be  executed  and  2)  from  Eqs  (4- 
If)  and  (4-18),  a  negative  value  of  kf  must  exist  to 

insure  a  right  turn.  Therefore,  kr  shall  be  selected 

and  kj  calculated  from  the  relationship 

- 

To  insure  that  a  negative  radical  is  not  chosen  ky  is 
limited  by  solving  Eq  (4-44)  for  ky  and  substituting 
k+  *0.  The  result  is 

'Gt£«b|s£l|  s  £  Gt£9"(s^r|  <4-46> 

It  is  of  interest  to  plot  Eq  (4-44)  to  aid  in  selecting  a 

reasonable  bank  angle  at  the  final  time.  Fig  4-1  is 
plotted  as  an  example  for  a  minimum  time-to-turn  problem 


Figure  4-1.  Minimum  Time  Final  Velocity  Contour  Plot 
for  selection  of  Ay  and  A^, 
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(Gt£-1)  for  final  velocities  of  500,  750  and  1000  fpa. 
Aa  seen  in  the  figure,  there  is  a  locus  of  maxima  for  all 
velocities  which  is  a  straight  line.  Solving  Eq  (4-44) 
for  the  value  of  Xy  which  maximizes  X*  yields. 


-  V 


g(N2-l) 


Afmax 


(4-47) 


The  values  of  Eq  (4-47)  in  turn  yield 


(4-48) 


which  is  the  tangent  of  the  bank  angle  for  a  level 
coordinated  turn  at  the  final  time.  Note  that  the  locus 
plotted  corresponds  to  82<>  of  bank  for  7.22g's  only, 
tf  Therefore,  choosing  combinations  of  Xy  and  X*  above 

the  locus  will  cause  the  aircraft  to  pitch  up  at  tf  and 
combinations  below  it  will  cause  a  pitch  down.  A  similar 
plot  can  be  generated  for  any  other  free  final  time 
problem. 

In  summary,  the  costate  final  conditions  are  stated 
as  follows x 

a)  minimum  time 

G  -  tf  A  -  [  0  0  0  0  XyX*  ]  (4-49) 

b)  maximum  energy 

G  -  Ei-Ef  A  »fo  0  -1  -V  xy  X*  1  (4-50) 

i  J 
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c.  maximum  energy  gradient 


G  .  Ei-Ef 


(4-51) 


Substitution  o f  table  4-1  quantities  into  Eq  (4-6)  yields 

Vc  ■  558.07yQI~  (4-52) 

From  Eq  (4-52),  it  is  clear  that  the  selection  of  a  final 
altitude  fixes  the  final  velocity.  Along  with  a  choice 
for  a  final  time,  the  values  of  Aft  and  Av  can  be 
determined  for  all  the  problem  formulations  investigated. 

For  the  selection  of  Ay  and  A*  for  fixed  final 
time  problems,  arbitrary  numbers  which  generate  a 
feasible  final  bank  angle'  are  the  best  guesses.  However, 
for  free  final  time  problems,  use  Eq  (4-52)  for  a 
velocity,  select  a  final  bank  angle  for  pitch  up  or  pitch 
down  at  tf,  pick  Ay  off  of  Fig  4-1  using  a  protractor 
and  substitute  it  into  Eq  (4-45)  to  calculate  A*  •  This 
process  will  always  insure  that  Hf  +  Gtf  ■  0.  The 
solutions  for  Ay  and  A#  for  both  the  high  and  low 
speed  minimum  time  cases  are  plotted  for  reference  on  Fig 
4-1. 

Once  these  final  conditions  at  tf  are  established, 
the  state  and  costate  differential  equations  (4-1)  and 
(4-34)  through  (4-39)  can  be  integrated  to  the  initial 
conditions  at  tQ.  The  methods  of  integration  and  search 
for  the  optimum  control  trajectory  are  covered  in  Sec  5. 
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V.  Numerical  Method  for  Solution  off  the  TPBVP 
As  stated  in  Section  4,  the  dependence  of  A^"  and 
Av“  on  Ay*  makes  forward  integration  of  the  state  and 
costate  equations  from  initial  to  final  conditions 
impossible.  Backward  integration  from  the  final  initial 
conditions  must  be  employed  for  both  minimum  time  and 
maximum  energy  problems. 

Minimum  Time-to-Turn  TPBVP  Solution 

Given  the  2n+l  dimensional  final  condition  vector* 
Zf,  having  the  quantities 


Zf  ■  [x  Y  h  V  y  ip  Ax  Ay  Ah  Xv  Ay  A*  H 

-  f - -  -  0  1800  o  0  0  0  v 2  -1  ]  (5-1) 


and  the  2n+l  dimensional  initial  condition  vector*  Zj, 
having  the  quantities 

Zi  «  X  Y  h  V  y  \p  Ax  Ay  Ah  Av  Av  Ay  H 

-  [  °  0  hi  Vi  0  0  -------1  (5-2) 


two  subset  vectors  relating  final  condition  variables  and 
initial  condition  requirements  can  be  generated. 

It  can  be  seen  from  Eqs  (4-1)  and  (4-34)  through  (4- 
39)  and  (4-41)  that  the  down  range  and  cross  range 
variables  X  and  Y  and  their  respective  influence 
functions  Ax  and  Ay  have  no  effect  upon  the  problem 
solution  as  they  are  dependent  variables.  They  can  be 
eliminated  from  the  solution  for  the  sake  of  expediency. 
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It  was  learned  from  Sec  4  that  a  selection  of  the 
altitude  at  tf  dictated  the  final  corner  velocity. 
Therefore#  velocity  at  tf  is  no  longer  a  variable. 
Additionally#  the  use  of  this  altitude- velocity  pair  in 
concert  with  the  requirement  for  Hf  +  Gtf  ■  0  allowed  for 
the  calculation  of  A*  at  tf  as  a  function  of  Ay  . 
Now#  A*  and  Hf  are  no  longer  variables.  As  this  is  a 
free  final  time  problem  with  an  interior  time  constraint# 
tf  and  tc  must  be  added  as  variables  to  be  selected 
before  the  integration  is  performed.  The  final  condition 
variable  subset  vector#  z#  which  will  be  varied  to  reach 
the  initial  conditions  is  therefore 

z  ■  f  h  Ay  tf  tc  j  (5-3) 

These  four  final  time  variables  match  with  the  four 
initial  time  requirements 

zj  -  [  hi  Vi  Vi  ]  (5-4) 


After  integrating  from  the  non-optimal  z  vector 
conditions  at  t«tf  to  t»t0,  there  will  be  a  non-zero 
residual  miss  vector  R  formed  by  the  difference  between 
the  initial  time  results  and  the  desired  conditions.  R 
is  generated  from  Z}  and  is 


R  -  [  h-hi,  V-Vi,  y  ,  \\)  ]  t-tc 


(5-5) 
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,sv>  The  objective  of  the  TPBVP  is  to  select  the  variables  in 

V* 

Z  which  when  integrated  from  tf  to  t0  drive  the  residual 
vector  R  to  zero.  In  order  to  reduce  the  values  in  R 
from  one  iteration  to  the  next#  the  direction  and 
magnitude  of  the  changes  to  the  Z  variable  values  must  be 
determined.  This  can  be  done  by  solving  a  Taylor  series 
expansion  of  R  about  the  z  vector.  From  (Ref  14:686)# 
this  is  found  to  be 

R(z+  6z)  ■  R(z)  +  dR  <5z  +  i  0zT  d2R  qz 

0  z  2  qz2 

+  Higher  Order  Terms  (5-6) 

Assuming  that  the  second  and  higher  order  terms  in  Eg  (5- 
6)  are  negligible  with  respect  to  the  first  order  term's 
influence  on  a  solution#  then  the  equation  for  change  in 
R#  A  R#  due  to  a  variation  in  the  final  time  vector  z# 
6  z,  is 

A  R  »  R(z+  6z)  -  R(z)  -  6z  (5-7) 

d  z 

It  is  desirable  that  after  a  single  variation  of  z, 
due  to  an  increment  6  z,  the  term  R(z+  <5z)  will  diminish 
to  zero.  In  order  for  a  vector  6z  to  be  calculated 
which  will  cause  R(z+  <5z)  to  be  zero#  Eg  (5-7)  is 
modified  to  be 

A  R  *  -  R(z)  ■  dR  6  z  (5-8) 

dz 
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Using  Newton's  Method  of  Tangents  (Ref  6:170)  to  solve  Eq 
(5-8)  for  6  z,  it  is  found  that 

<$2  -  -  3R"1  R(z)  (5-9) 

d  z 

P.(z)  is  a  straight-forward  numerical  difference  as  shown 
in  Eq  (5-5).  The  4x4  first  order  derivative  matrix  is 
formed  using  a  numerical  central  difference  scheme  (Ref 
6s220)  operating  on  each  of  the  four  variables  in  z. 
This  is  then  inverted  and  used  to  calculate  6  z  in  Eq 
(5-9). 

6 z  is  the  vector  along  which  a  search  is  conducted 
to  find  what  fraction  of  its  magnitude  will  cause  the 
biggest  reduction  is  the  R(z)  vector.  Should  a  full 
magnitude  search  along  6z  not  reduce  the  R  vector  on 
the  first  attempt,  reductions  in  the  magnitude  of  <5  z 
and  subsequent  searches  are  required  until  the  maximum 
reduction  in  R  is  reached.  The  equation  for  reduction  in 
6  z  is  formed  using  the  Golden  Section  Ration  (Ref 
6:460)  and  is 

0zj+l  -  0.618034  •  dzj  j>l  (5-10) 

Reductions  in  <5z  are  required  until  the  maximum  reduction 
in  R  is  achieved. 

In  summary,  the  following  steps  are  taken  tc  obtain 
the  optimal  solution  for  the  minimum  time  problem: 
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1)  Select  values  for  the  elements  of  z  which  comply 
with  the  transversality  conditions  of  Section  4. 

2)  integrate  the  state  and  costate  equations 
backwards  constrained  by  the  requirements  of 
Section  4  from  tf  through  tc  to  t0  and  form  the 
residual  vector  R. 

3)  Numerically  form  the  derivative  matrix  JLR  by 

d  z 

integrating  backward  to  find  the  variations  in  R 
due  to  independent  variations  in  the  z  elements. 

4)  Calculate  6 z  from  Eqs  (5-9). 

5)  Modify  the  values  of  z  by  6z  and  perform  a  search. 

6)  If  a  reduction  in  R  does  not  occur  on  this 
search,  decrease  6z  according  to  Eq  (5-10)  and 
repeat  step  5  until  a  reduction  occurs. 

7)  If  the  root  sum  square  (RSS)  of  elements  of  the 
reduced  R  vector  is  greater  than  10-4,  then 
repeat  steps  2  through  7. 

8)  If  the  RSS  of  R  is  less  than  10"4,  then  the 
solution  has  been  achieved  for  the  minimum  time 
problem. 


Maximum  Energy  Turn  TPBVP  Solution 

The  maximum  energy  turn  formulation  is  nearly 
identical  to  the  minimum  time  formulation.  As  Eq  (4-43) 
does  not  apply  to  a  fixed  final  time  problem  such  as 


this,  there  are  only  2n  final  conditions  in  the  Z  vector. 
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The  resultant  modifications  to  the  a  vector  are  that  tf 
is  no  longer  a  variable  and  X*  has  taken  its  place. 
This  still  leaves  4  variables  and  4  equations  to  solve 
and  the  same  methods  used  for  solving  for  da  can  be 
used. 

Maximum  Energy  Gradient  TPBVP  Solution 

The  method  used  to  solve  this  formulation  is 
identical  to  the  minimum  time-to-turn  method.  Only  the 
constant  value  for  the  Hamiltonian  has  changed  due  to  the 
new  values  for 
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VI.  Anticipated  Results 

No  predictions  for  the  values  of  tirade*  Etrade  or 
EGmax  can  be  made  other  than  the  qualitative  predictions 
on  Figures  2-1  and  2-2.  However,  qualitative  predictions 
about  the  control  histories  can  be  presented.  They  are 
presented  in  Pigures  6-1  through  6-6. 

Optimal  Control  Histories 

The  predicted  bank  angle  histories  for  Vf  >  Vc  are 
shown  in  Pig  (6-1).  Por  the  minimum  time  turn  (curve 
'a'),  an  initial  bank  angle  of  less  than  90°  is  required 
to  climb  while  turning  to  bleed  off  airspeed  as  rapidly 
as  possible.  Subsequent  to  reaching  the  corner  velocity 
arc.  the  bank  angle  is  increased  to  be  greater  than  90° 
to  restore  the  flight  path  angle  end  condition.  Given 
more  time  in  which  to  complete  the  turn,  a  slight 
decrease  in  the  initial  bank  angle  and  a  slight  increase 
in  final  bank  angle  should  be  in  order  as  climbing  for  a 
longer  duration  to  a  higher  corner  velocity  should  occur 
with  resultant  greater  flight  path  angle  excursions. 
However,  the  basic  shape  of  the  history  should  remain 
intact. 

The  lift  coefficient  histories  shown  in  Fig  (6-3)  all 
represent  maximum  load  factor  flight  until  the  corner 
velocity  arc  is  encountered  at  which  point  CLmax  is  the 
optimum  control.  An  interior  lift  coefficient  may  occur 
should  time  be  available,  and  insufficient  climb  and 
acceleration  capabilities  exist  at  the  higher  altitudes 
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at  the  maximum  load  factor. 

The  throttle  history  for  Vi  >  Vc  flight  is  the  most 

provocative  control  to  discuss  for  this  problem.  Prom 
(Ref  1)#  it  is  known  that  the  throttle  setting  will  be 
minimum  until  the  corner  velocity  arc  is  intercepted  at 
which  point#  an  interior  thrust  control  will  be  employed 
(curve  'a')*  Given  extra  time  to  complete  the  turn  and 
maximize  the  final  energy  it  is  postulated  that  the 
throttle  will  initially  be  at  minimum  setting  and  then 
switch  for  some  time  to  maximum  throttle  prior  to  an 
interior  corner  velocity  arc  setting.  This  would  be  a 

" bang-bang-singular"  throttle  control  (curve  'b' ).  This 

combination  is  proposed  as  it  would  require  less  power 
and  time  to  decelerate  from  high  speed  and  then 
reaccelerate  to  corner  velocity  arc  than  the  converse. 
Given  enough  time#  no  minimum  control  bound  would  be 
expected;  only  maximum-singular  control  (curve  'c')  as 
will  be  discussed  for  >VC. 

For  the  case  where  Vj  <  Vc,  the  control  histories 
should  differ  quite  markedly  from  the  V£  >  Vc  case.  In 
Fig  (6-2),  the  initial  bank  angle  for  the  minimum  time 
turn,  curve  'a',  should  be  the  largest  of  all  the  turns 
so  as  to  create  the  greatest  acceleration  possible 
terminating  in  the  smallest  bank  angle.  Given  more  time 
to  turn,  (curve  '  b’  ),  less  diving  and  more  climbing  would 
result,  thus  reducing  the  bank  angle  requirements.  If 
sufficient  time  and  excess  power  are  available,  the  bank 
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angle  may  start  at  a  value  less  than  90°  and  finish 
greater  than  90°  to  maximise  final  altitude  and  airspeed. 

The  lift  coefficient  histories,  as  seen  in  Fig  (6-4), 
should  all  be  at  maximum  lift  coefficient  unless 
sufficient  time  to  complete  the  turn  is  supplied  wherein 
a  corner  velocity  arc  is  not  part  of  the  optimal  history* 
In  that  case,  an  interior  lift  coefficient  may  occur,  but 
is  not  expected  in  this  exercise. 

From  Fig  (6-6),  the  throttle  histories  for  Vj <  Vc  all 
contain  maximum  throttle  settings  until  the  corner 
velocity  arc  is  reached  (curve  'a'),  otherwise,  known  as 
"bang-singular"  control.  As  time  to  turn  is  increased, 
climbing  and  accelerating  to  higher  speeds  should  cause  a 
larger  tc  value  and  a  higher  corner  arc  throttle  setting 
(curve  'b').  Given  sufficient  time  to  turn,  a  corner 
velocity  arc  may  not  be  part  of  the  optimal  history  and 
the  result  is  a  maximum  throttle  setting  throughout  the 
turn.  The  shape  of  the  throttle  history  during  corner 
velocity  flight  is  dictated  by  an  altitude  rate  and 
flight  path  angle  rate  combination. 

None  of  the  plots  presented  are  intended  to  be 
construed  as  predictions  of  exact  values  in  either  time 
or  control  variable.  They  are  intended  to  predict  trends 
for  comparisons  between  minimum  time  and  maximum  energy 
turns. 
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Figure  6-3.  Lift  Coefficient  vs.  Time  for  V. >V 
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Figure  6-4.  Lift  Coefficient  vs.  Time  for  V.<V 
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Figure  6-5.  Throttle  Setting  vs.  Time  for  V. >V 
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Figure  6-6.  Throttle  Setting  vs.  Time  for  V.  <V 
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VII.  RESULTS  OBTAINED 

The  equations,  conditions  and  constraints  for  all 
three  problem  formulations,  set  forth  in  sections  4  and 
5,  were  programmed  on  a  32  bit  DEC  20/20  digital  computer 
in  the  Space  Division  Computer  Center.  A  copy  of  the 
computer  programs,  OPTMYZ  and  its  subroutines,  are 
located  in  the  appendix.  The  results  of  this  program's 
execution  are  presented  in  this  section. 

Minimum  Time-to-Turn  Results 

As  the  costate  variables  determine  the  optimal 
controls  which  in  turn,  determine  the  state  trajectories, 
the  results  are  presented  in  that  order.  Time  history 
plots  for  the  costates,  controls  and  states  are  presented 
in  figures  7-1  through  7-9.  By  examining  these  plots, 
several  observations  can  be  made. 

First,  the  costate  histories  figs  (7-1)  through  (7-3) 
are  well  behaved  and  are  primarily  monotonic;  the 
discontinuities  in  Ah  and  Av  do  not  create  erratic 
behavior  or  unexpected  characteristics.  It  may  be 
coincidence  and  of  no  consequence  but  it  is  of  interest 
to  note  that  for  Ah  when  t  <  tc,  is  virtually  a  constant 
while  is  on  the  Nmax  boundary  for  the  high  speed 

solution  but  varies  smoothly  when  C l  is  equal  to  CLmax 
for  the  low  speed  solution.  Otherwise,  Ah  acted  as 
expected.  Av  responded  per  the  optimal  control 

requirements  of  Eq  (4-24)  through  Eq  (4-32).  It 
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Figure  7-1.  A  h  vs  Time  for  Minimum  Time  High 
and  Low  Speed  Solutions 


Figure  7-2.  Av  vs  Time  for  Minimum  Time  High 
and  Low  Speed  Solutions 
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exhibited  no  unexpected  trends.  Figure  7-3  shows  the 
correlation  between  the  bank  angle,  flight  path  angle, 
Ar  and  the  constant  A*  .  The  bank  angle  history  can  be 
determined  from  the  plot  using  a  protractor  the  same  way 
as  was  done  in  fig  (4-1)  and  can  be  correlated  roughly  to 
time  using  the  tick  marks  at  one  second  increments. 
Again,  these  parameters  all  vary  smoothly  with  time. 

The  control  histories  are  shown  in  figures  (7-4) 
through  (7-6).  The  bank  angle  history  shapes  were  as 
expected  although  the  amplitudes  were  larger  than 

anticipated.  The  lift  coefficient  histories  were  exactly 
as  expected  for  both  high  and  low  speed  solutions.  The 

shapes  of  the  throttle  histories  were  as  expected, 
although  like  the  bank  angle,  their  excursions  were 
greater  than  anticipated.  As  seen  in  fig  (4-6),  the 

singular  arc  portion  of  the  solutions  appear  to  go  right 
up  to  the  boundary  denoted  by  the  arrow  heads.  In 
actuality,  the  values  of  jt  at  these  two  points  are  n  • 
0.9999  at  t  ■  6.12  secs  for  the  high  speed  solution  and 
n  »  -0.4997  at  t  »  6.22  sec  for  the  low  speed  solution 
thus  keeping  them  interior  to  the  boundaries.  However, 
in  order  to  insure  that  it  remained  interior,  a  special 
search  routine  in  the  integrator  had  to  be  programmed  to 
look  for  the  point,  if  one  existed,  where  the  arc  touched 
the  boundary,  thus  ending  the  arc.  More  will  be  said 
about  this  in  the  technique  problems  section.  In 

retrospect,  it  makes  sense  that  the  throttle  would  go 


Bank  Angle,  0  (DEG) 
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Figure  7-4.  Bank  Angle  vs  Time  for  Minimum  Time 
High  and  Low  Speed  Solutions 
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Figure  7-5.  Lift  Coefficient  vs  Time  for  Minimum  Time 
High  and  Low  Speed  Solutions 
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Figure  7-6.  Throttle  Setting  vs  Time  for  Minimum  Time 
High  and  Low  Speed  Solutions 
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right  up  to  the  boundary  oven  on  tha  aro  bacauaa  the 
control*  will  use  all  of  tha  authority  thay  hava  to 
minimize  tha  tima  to  turn.  Just  a*  in  a  bang-bang 
controller,  a  bang-singular  controller  should  ba  axpactad 
to  go  from  boundary  to  boundary  while  meeting  its 
interior  control  constraints. 

Turning  to  the  state  trajectories  starting  with  Fig 
(7-7)  it  is  seen  that  both  solutions  end  up  on  the 
singular  arc  and  both  met  a  small  loss  in  energy  (-262  ft 
for  low  speed  and  -1258  for  high  speed).  Both  maneuvers 
are  basically  planar,  the  steepness  of  which  are  dictated 
primarily  by  the  initial  bank  angle  values.  This 
steepness  correlates  directly  to  the  flight  path  angle 
histories  of  fig  (7-8),  the  low  speed  being  more  vertical 
a  maneuver  than  the  high  speed.  The  main  reason  for  them 
not  both  being  vertical  plane  maneuvers  is  that  the 
throttle  settings  on  the  arc  were  already  up  against  the 
boundary.  However,  increasing  the  thrust-to-weight  ratio 
to  approximately  1.9  or  greater  would  result  in  two 
vertical  maneuvers  where  the  arc  throttle  histories  would 
not  be  nearly  so  close  to  the  boundaries  because  the  use 
of  gravity  would  be  maximized.  The  heading  angle 
histories  show  nothing  unexpected.  The  low  speed 
solution  has  higher  peak  turn  rates  due  to  the  steeper 
flight  path  angle  history.  Tables  (7-1)  and  (7-2) 
summarize  the  key  values  for  the  states  and  costates  for 
the  low  and  high  speed  solutions  respectively.  From  both 
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Figure  7-8.  Flight  Path  Angle  vs  Time  for  Minimum 
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Table  7-1.  End  and  Corner  Values  of  States  and  Costates 
for  Low  Speed  Minimum  Time  Solution 
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t^-10.87  t  -4.  495-t,. 

f  c  c 


*0-° 


0.0002078  0.0001103  0.0001101 


0.008138 


0.005010 


2.120 


-1.917 


3.108 


-2.1075 


-2.1075 


-2.1075 


1239. 5 


3252.2 


3082.0 
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■1251.8 


-4668.3 


19033 


16372 


14998.3  rh--1.7438 


762.6 


730.0 


960.05 


ry-0.3734 


0.7714 


0.006517  r  “0.006517 


3.14159 


0.8626 


0.06289  r  =0.06289 


||R|| -1.7555 


Table  7-2.  End  and  Corner  Values  of  states  and  costates 
for  High  Speed  Minimum  Time  Solution 
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tables,  it  is  seen  that  the  end  condition*  were  not  net 
to  the  accuracy  that  was  expected*  The  norma  of  the 
residual  vectors  shown  are  as  low  as  could  possibly  be 
achieved  with  the  program,  computer  and  methodology 
employed. 

Problems  with  the  Technique 

The  search  for  the  optimal  controls  began  by  using  a 
variable  stepsise  integrator  and  allowing  only  small 
increases  in  the  total  time  allotted  for  the  turn  to 
occur  on  each  iteration  to  insure  stability  in 
approaching  the  solution.  However,  after  achieving  a 
norm  of  the  residual  vector  of  approximately  10,  fixed 
steps! ze  integration  had  to  be  employed.  The  prime 
reason  for  this  is  due  to  the  special  programming  which 
had  to  be  added  to  the  integrator  to  find  the  point  in 
time  where  7r  went  so  close  to  or  touched  the  boundary 
during  the  search  phase.  This  search  could  not  be 
performed  using  variable  stepsize  integration. 

The  use  of  small,  fived  stepsize  integration  helped 
to  get  closer  to  the  optimal  solution,  but  the  extreme 
sensitivity  of  the  norm  of  the  residual  vector  to  the 
final  value  of  Ay  created  havoc  with  finding  the  minimum 
in  the  state  space.  As  shown  in  fig  (7-10),  it  is  seen 
that  the  standard  parabolic  minimum  seen  in  bang-bang 
minimum  time  problems  has  been  replaced  with  a  gross 
discontinuity  in  the  norm  at  the  point  where  the 
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Figure  7-10.  Norm  of  Residual  Vector  vs  Ay  for 

Minimum  Time  High  and  Low  Speed  Solutions 
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magnitude  of  Ay  gets  large  enough  to  create  a  steeply 
banked  planar  maneuver  which  causes  ir  to  require  a 
boundary  value  at  the  arrow  head  points  on  fig  (7-6). 
This  causes  the  time  at  which  the  switch  from  the 
singular  arc  to  the  boundary  value  control  to  snap  back 
and  forth  during  the  search  for  the  minimum.  In  the  low 
speed  case#  the  change  in  Xy  from  -3.677  to  -3.678 
caused  a  discontinuous  jump  in  the  norm  from  less  than 
1.0  to  more  than  100.  whether  the  sensitivity  matrix  for 
calculating  6z  was  formed  using  a  forward#  central  or 
backward  difference  scheme  was  immaterial.  The  z 
components  being  so  tightly  coupled  could  not  be  modified 
so  slightly  in  the  proper  bz  direction  to  reduce  the 
norm  any  further  than  that  shown  in  tables  (7-1)  and  (7- 
2). 

Even  though  the  end  conditions  were  not  met  close 
enough  to  form  a  residual  norm  of  10“4  or  less#  the 
trajectories  shown  are  extremely  close  to  the  optimum  and 
demonstrate  that  the  technique  and  problem  formulation 
are  correct  and  feasible. 

Maximum  Energy  Turn  Results 

No  feasible  trajectories  were  generated  for  the 
maximum  energy  turn  formulation  which  came  close  to  the 
desired  end  conditions.  No  norm  of  the  end  conditions 
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single  set  of  final  time  conditions  which  would  generate 
one  feasible  trajectory  from  which  to  propagate  the 

others.  Several  fixed  final  times  were  used  as  starting 
points,  but  all  times  generated  approximately  the  same  end 
results.  The  multipliers  Xy  and  X*  grew  an  order  of 
magnitude  with  each  iteration  surpassing  values  of  1013 
and  growing  with  no  reasonable  reduction  in  the  end 

condition  norm.  No  definite  reason  can  be  given  for  this 
type  of  performance;  perhaps  the  linearity  of  the 

aircraft  thrust  and  drag  characteristics  create  a  flat 
state  space  where  there  is  no  obvious  maximum  energy 
turn. 

Since  no  maximum  energy  turn  solutions  were  achieved, 
the  graphical  method  of  fig  (2-1)  for  depicting  the 

energy-time  tradeoff  can  not  be  used. 

Maximum  Energy  Gradient  Results 

Unfortunately,  the  same  type  of  results,  as  in  the 
maximum  energy  turn  attempt,  occurred  here  as  well.  The 
values  for  Xy  and  X*  were  in  the  vicinity  of  108  and 
growing  steadily  showing  no  reduction  in  the  norm  of  the 
end  conditions.  Additionally,  the  final  time  would  not 
settle  in  one  neighborhood  but  would  grow  and  decrease  as 
the  attempts  progressed.  Once  again,  no  positive 
explanation  can  be  offered  other  than  presented  in  the 
previous  section. 

No  maximum  energy  gradient  solution  was  achieved 
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using  this  method  either.  All  three  problems  used  pure 
optimal  control  formulations  and  ran  into  extreme 
difficulty  in  achieving  a  solution.  Evidently,  these 
formulations  do  not  lend  themselves  well  to  methods  using 
Lagrange  multiplier  techniques  and  may  be  more  well 
suited  for  parameter  optimization  formulations. 
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VIII.  CONCLUSION 

The  goal  o£  this  effort  was  to  define  a  tradeoff 
between  time-to-turn  and  energy  in  a  turn  using  the 
classical  optimal  control  formulation.  First,  a  minimum 
time-to-turn  trajectory  was  to  be  generated  for  initial 
conditions  of  equal  energy  both  above  and  below  the 
corner  velocity  for  that  energy.  Second,  maximum  energy 
trajectories  beginning  at  those  same  initial  conditions 
were  to  be  generated  for  final  times  greater  than  the 
minimum  time.  From  all  of  these  trajectories,  the 
maximum  energy  gradients,  defining  the  time-energy 
tradeoffs,  were  to  be  graphically  identified.  Thirdly, 
these  tradeoffs  were  to  be  corroborated  using  the  maximum 
energy  gradient  formulation. 

This  project  was  successful  in  defining  the 
equations,  conditions  and  constraints  required  to  reach 
the  above  goals.  The  computer  program  OPTMYZ  was  written 
and  executed  on  a  DEC  20/20  digital  computer  which  was 
successful  in  defining  the  two  minimum  time  trajectories 
for  the  initial  conditions  above  and  below  the  corner 
velocity  conditions.  The  results  are  presented  in 
section  7. 

This  effort  failed  to  produce  feasible  maximum  energy 
turns  or  the  maximum  energy  gradients  required  for  the 
tradeoff  definition.  Therefore  the  time-energy  trade  can 
not  be  identified  for  either  of  the  initial  conditions 
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for  the  aircraft  model  used  and  the  formulations 
employed. 

Recommendations  for  further  work  in  this  area  are 
presented  in  section  9. 
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IX.  RECOMMENDATIONS 

The  following  are  prioritized  recommendations  for 
further  research  in  the  area  of  time-energy  tradeoffs. 

1)  Continued  effort  should  be  placed  upon 

defining  the  time-energy  tradeoffs  as  defined 
in  section  2  using  the  formulations  presented 
herein.  Emphasis  should  be  placed  upon 
developing  a  technique  which  will  better 
determine  the  minimum  cost  function  for  a 
complex  singular  arc  problem  such  as  the 
minimum  time  problem  presented. 

2)  The  time-energy  tradeoff  should  also  be 

determined  as  in  recommendation  (1)  using  a 
parameter  optimization  formulation  to  study 
the  tradeoff  of  complexity  versus  accuracy. 
Based  upon  these  results,  the  remaining 
recommendations  would  use  the  superior 
technique  for  future  work. 

3)  Recommendations  (1)  and  (2)  should  then  be 

repeated  using  a  more  realistic  aircraft 
model.  There  should  remain  regions  in  the 
altitude  velocity  map  where  high  specific 
power  exists,  but  the  thrust  should  vary  as  a 
function  of  altitude  and  Mach  or  Reynolds 
number  as  should  the  lift  and  drag 

parameters.  cI*max'  CD0  and  k*  A  thrust 
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pinch,  dynamic  pressure,  Mach  number  and 
altitude  placards  could  be  added  for  realism. 

4)  Recommendation  (3)  should  be  repeated  with 
the  addition  of  a  fuel  flow  or  varying  weight 
and  where 


w  ■  -c|tt| 


(9-2) 


5)  Recommendation  (4)  should  be  repeated  with 
the  performance  index  modified  to  pay  a 
penalty  for  excessive  fuel  consumption.  That 
is 


G  x(tf),tf 


hf-h!  +  V£2-V!2 

(Wf-Wjltf  2g(Wf-Wi)tf 


This  effort  would  be  the  most  direct  way  of 
looking  at  trading  time  and  fuel  for  energy 
with  the  most  realistic  model  for  an 
aircraft. 

6)  Using  the  results  from  the  previous 
recommendations,  a  quasi-iinear  feedback  law 
could  be  developed  for  two  modes  of 
operation.  The  first  mode  would  be  an 
algorithm  for  minimum  time  turn  guidance  free 
from  the  use  of  influence  functions  and  the 
second  mode  would  be  for  a  maximum  time- 
energy-fuel  tradeoff  based  upon  only  state 
and  state  rate  feedback. 
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7)  The  culmination  of  this  series  of  efforts 
would  be  to  program  these  feedback  laws  into 
a  manned  simulator  with  visuals  and  head-up 
displays  to  examine  the  man-machine 
applicabili ty. 


S  - 
|uV.\ 
k,V.\ 


.  X-' 


64 


lAiHuk 


AFIT/GAE/AA/78D-5 


Bibliography 


1.  Anderson,  Gerald  M.  and  W. L.  Othling  Jr.  "Optimal 
Trajectories  of  High  Thrust  Aircraft,"  Journal  of 
Aircraft,  13 x  180-184  (March  1976). 

2.  Bryson,  A. E. ,  W.F.  Denham  and  S.E.  Dreyfuss. 

"Optimal  Programming  Problems  with  Inequality 
Constraints,"  AIAA  Journal,  It  2544-2550  (Nov 

1963).  ~ “ " 

3.  Bryson,  A.E.  and  Y.C.  Ho.  Applied  Optimal  Control 

(Second  Edition).  Washington,  D.C. :  Hemisphere 

Publishing  Company,  1975. 

4.  Hedrick,  J.K.  and  A.E.  Bryson.  "Three  Dimensional 
Minimum  Time  Turns  for  a  Supersonic  Aircraft," 
Journal  of  Aircraft,  9t  115-121  (February  1972). 

5.  Humphreys,  Robert  P.  ,  G.R.  Hennig,  W.A.  Bolding  and 

L.A.  Helgeson.  "Optimal  Three  Dimensional  Minimum 
Time  Turns  for  an  Aircraft,"  Journal  of 

Astronautical  Sciences,  20:  88-112  (September- 

October  1972). 

6.  Ketter,  R.L.  and  S.P.  Prawel  Jr.  Modern  Methods  of 
Engineering  Computation.  New  Yorks  McGraw-Hill  Book 
Company,  196$. 

\ 

7.  Miele,  Angelo  A.  "Recent  Advances  in  Gradient 
Algorithms  for  Optimal  Control  Problems,"  Journal  of 
Optimization  Theory  and  Applications,  17 t  361-430 
(December  1975). 

8.  National  Oceanic  and  Atmospheric  Administration. 
U. S.  Standard  Atmosphere,  1976.  Washington,  D.C. : 
National  Aeronautics  and  Space  Administration,  United 
States  Air  Force,  1976. 

9.  Newman,  William  M.  and  R.F.  Sproul.  Principals  of 
Interactive  Computer  Graphics.  New  York:  McGraw-Hill 
Book  Company,  19*73 . 

10.  Nicolai,  Leland  M.  Fundamentals  of  Aircraft  Design. 
Fairborn:  E.P.  Domicone  Printing  Services,  197!>. 

11.  Noble,  Ben  and  J.W.  Daniel.  Applied  Linear  Algebra. 
Engelwood  Cliffs:  Prentice-Hall,  Incorporated,  1977. 

12.  Rader,  James  E.  Lecture  Notes  distributed  during  the 

Optimization  Techniques  Sequence.  School  of 

Engineering,  Air  Force  Institute  of  Technology, 
Wright-Patterson  AFB,  1978. 


AFIT/GAB/AA/78D-5 


13.  Rader,  James  B.  and  D.G.  Hull.  "Computation  of 

Optimal  Aircraft  Trajectories  using  Parameter 
Optimization  Methods,"  Journal  of  Aircraft.  1 2 » 
864-866  (November  1975).  "" 

14.  Wylie,  C.R.  Advanced  Engineering  Mathematics.  New 
York*  McGraw-Hill  Book  Company,  1966. 

15.  Yajnik,  K.S.  "Energy-Turn  Rate  Characteristics  and 

Turn  Performance  of  an  Aircraft,"  Journal  of 
Aircraft,  14*  428-432  (May  1977). 


AFIT/GAE/AA/78D-5 


Appendix 

Computer  Program: 
OPTMYZ 
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program  optmvz 

dimension  x(18).xf(19).xi(19) ,dx(4),r(4),rm(4),rp(4),dz(4), 
&drdz(4,4),drdzi(4,4),d(4,4) 
real  k,nmax 

common  a,beta,olmax,edo,g,k,nmax,tmax,w,xi4 
oommon  olnmax ,hfifi, Jflag,qa, rosig 
data  xi<6),xi(7>,xf(6),xf(7)/«.8,8.9,8.8,3.14159/ 
data  a, beta,olmax,cdo, g,k,nmax,tmax ,  w/ 

4237. 9, «9. 88983281 1628, 1.9,8, 92, 32. 174,9. 95 ,7.22, 

&18225.9, 12159.9/ 
write(5 , 1 ) 

1  format(1x, 'load  lflag  for  tmin,Emax  or  Edotmax') 
read(5 , *)  lflag 

write(5,2) 

2  format( 1 x ,’ load  hf ,lambdagammaf ,lambdapsif , tf , tc ,hl , vi ' ) 
read (5,*)  xf (1 ) ,xf (2) ,xf(3) ,tf ,tc,xi(1 ) , x i ( 4 ) 
xi4=xi(4) 

write<5,3) 

3  format ( 1 x do  you  wish  to  print  each  step(yes|9,no»1 )?' ) 
read(5,*)  iprint 

call  assign(xf,xi,tf, lflag) 
do  4  i  =  1 , 19 

4  x(i)=xf(i) 

call  integ(19,x,tf,tc, iprint) 

xf(9)=9x(9) 

xf ( 18)=«x( 18) 

call  ndmiss( x , xi ,r ,rnorm) 

write(5,5)  r,rnorm 

5  format ( 1 x , ' r  and  rnorm ’ , / ,5 ( 1 x , f 19.4 ) ) 
it=9 

6  it=it+1 
write(5,7) 

7  format ( 1 x do  you  want  to  continue?') 
read(5,*)  choice 

if ( choice, eq.8 )  go  to  16 
if(rnorm,le, 1.9898)  go  to  16 
do  12  i= 1 , 4 
i  i  =  i 
iii=1 

do  8  j  =  1 ,19 

8  x(j)sxf(j) 
tfsstf 
tcs=tc 
write(5,81 ) 

81  format(1x, 'deltas') 

call  deltas( x ,dx , tfs , tes , lflag , ii , iii) 
write(5,81 ) 

„  write(5.82) 

82  formatUx,  'assign') 

call  assign(x,xl, tfs, lflag) 
write(5, 82) 
write(5 , 83 ) 

83  format (1 x ,' integ ' ) 
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call  integ(1fl,x,tfs,tcs,iprint) 

write(5,83) 

write(5,84) 

84  format(1x, 'ndmiss' ) 

call  ndmiss(x,xi,rp,rnormp) 
write(5, 84) 
write(5,9)  rp,rnormp 

9  format(1x,  'rp  and  rnormp' ,/,5(1x,f18.4)) 
iiiMl 
do  18 

19  x(j)=xf(j) 
tfsatf 
teas tc 

call  deltas ( x ,dx , tfs , tes , If lag ,  ii , ill ) 
call  assign(x,xi, tfs,lflag) 
call  integ(18,x,tfs,tcs,iprint) 
call  ndmiss(x, xi , rm, rnormm) 
write(5,11)  rm, rnormm 

11  format(1x, 'rm  and  rnormm' ,/,5( lx ,f 18,4) ) 
do  12  js1,4 

12  drdz(j,i)=(«‘p(J)*rm(j))/(2.9),dx(i)) 
write(5,13)  dx 

13  format ( 1 x, 'dx',/,4(1x,f14*7),/) 
call  mtxinv( drdz , drdzi ,d ,4 ) 

do  14  i=  1 , 4 
dz(i)=8.8 
do  14  j  =  1,4 

14  dz(i)sdz(i)|drdzi(i,j)#r(j) 
write(5,15)  dz 

15  format (lx, 'dz' ,/,4(1x,f18#4),/) 

call  search(xf , xi , tf , tc ,dz , r , rnorm, lflag , ipr int) 
write(5,17)  xf,tf,tc,r, rnorm 
if (rnorm.gt, 1 .9998)  go  to  6 

16  write(5,17)  xf ,tf ,tc ,r , rnorm 

17  format(1x, ' xf ,tf , tc ,r , rnorm' ,/,2(5(1x,f13.7),/), 
&2(1x,f19.6),/,5(1x,f13.7>) 

if(choice,eq.9)  go  to  28 
do  18  i*1 , 18 

18  x(i)sxf(i) 

call  integ( 18 , x , tf , tc  ,9 ) 
call  ndmiss(x,xi,r, rnorm) 
wrlte(5,19)  x,tf ,tc ,r , rnorm 


29  stop 
end 
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subroutine  assign( x , xi , tf ,lflag) 
dimension  x(18),xi(l8) 
rsdl  nmax 

data  beta,g,nmax/§9. 88883281 1628,32. 174,7.22/ 

x(4)=558.87441*sqrt(exp(«beta*x(1))) 

if (lflag.eq.2)  go  to  18 

Eisxi(1 )+xi(4)«*2/2.8/g 

Esx(1)+x(4)*#2/2.8/g 

if ( lflag.eq. 1 )  Gtf=1.8 

if ( lflag.eq.3)  Gtf=(E#Ei)/tf»*2 

x2min=*Gtf*x(4)/288.123 

x2max~Gtf*x(4)/264.471 

if(x(2).lt.x2min)  x(2)=x2min 

if ( x(2) .gt.x2max )  x(2)=x2max 

x(3)=«sqrt((x(4)/g/nmax)»*2*(Gtf#x(2)*g/x(4))**2«x(2)*»2) 

x(5)=8.8 

x(8)s8.8 

if ( lflag.eq. 1 )  go  to  29 


19  x(5)=«1.8 
x(8)s«x(4)/g 
if (lflag.eq.2)  go  to  28 
x(5)=x(5)/tf 
x(8)=x(8)/tf 
28  return 
end 
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subroutine  integ(n,x,tf , to ,i print) 
dimension  x(1i) 
real  k,nmax 

common  a,beta,clmax,cdo,g,k,nmax,tmax ,w,xi4 

common  clnmax ,hfifi , jflag,qa,rosig 

Jflag=8 

dts|8 *25 

tcs=tc 

call  rk78(n,x,tf ,tcs,dt ,tc , iprint) 

x(5)=x(5)*x(8)*beta»x(4)/2.« 

x(8)=fl.i 

tcxstcs«1 «9et96 

dta««.25 

call  rk78(n,x,tcx ,fl.9,dt , tcs, i print) 

testes 

return 

end 


AFIT/GAB/AA/78D- 5 


subroutine  rk78(n,x,t1 ,t2,dt, tc,i print) 
dimension  alph( 13) ,b( 13, 12) , ch( 13) . x( 18 ) . xdum( II ) 
dimension  fit  18)  tf£( 18 )  ,f3( fa) , f*z*<  fa ) ,  f5(  1«)  ,f6(  If ) , f7(1 8 ) , 

&f8(18) ,f9(18) ,f18(18) ,f11(18),f12(18) ff13(1i) 
real  k,nmax ,mu2,mu3 

common  a, bet a ,clmax ,cdo,g,k ,nmax ,tmax ,w, xi4 
common  clnmax ,hfif i , jflag , qa , rosig 
data  (alph(i) ,i=1 , 13)/ 

&8.  ,8.748748748748748e*81  ,8.11 111111 1111 11 1e+88, 

&8. 166666666666666e+88 ,8.41 666 6 6 66666666 e+88 ,9*5  e+88, 

&8 .8333333333333336+88 , 8. 166666666666666e+88 ,8.666666666666666e+88 , 
&8 • 33333333333 3333e+88 ,1.8  e+88, 8. 

&1.8  e+88/ 

data  ( (b(i , j) , i=1 , 13) , j=1 ,3)/ 

&8.  ,8.74874874«748748ei81 ,i.277777777777777e«81 , 

&8.416666666666666ef81,8.416666666666666e+88,8.5  e«81 , 

&«. 231 481 481 481 481 e+88, 8. 183333333333333e+88 ,2. 8  e+8fl, 

&*.842592592592592e+88, 8*58121 95 121 95 121 e+88 ,8. 14 63 41463414634 ei 81, 


&8. 43 341 46341 46341 e+88 

,8. 

,8. 

&8.833333333333333ea81 

» 8  • 

?  8  • 

&8. 

*  8. 

,8. 

&8. 

» 8  * 

,8. 

&8. 

,8. 

,8. 

&8. 

,8. 

,8.125 

&«. 15625  e+81 

,8. 

>8. 

&8. 

,8. 

,8. 

&Q. 

*8. 

1 8  • 

data  ( (b( i , J) , i=1 , 13) 

» j=4,6)/ 

&8. 

» 8  . 

,8. 

e+88 


&Q.  ,8.15625  e+81,8.25  e+88, 
&8. 1 1574 874874 874 8 e+81 ,8.  ,§.883333333333333e+81 , 
&8, 21 2962962962963 e+88 , 8.287 92682926 8293 e+81 ,8.  , 
&«.287926829268293e+81 ,8,  ,8.  , 
&8 .  ,8.  ,8.  , 
&8.2  e+88, 8.248748748748741 e+81 ,8,271 1111 11 111111 e+«8 , 
&8.156444444444444e+82, 8. 722962962962963e+81, 8. 4386341 46341 463e+81, 
&Q.  , 8.4386341 46341 463 e+81 ,8.  , 
&8.  ,8.  ,8.  , 
&8.  ,8.  ,8.231481481481481 e+81 , 
&f. 222222222222222 e+88 ,8.11 8888888888888e+82 , 8.575925925925926e+81 , 
&8. 3678731 78731 787e+81 ,8. 146341 46341 4634e+88, 8. 352439824398244e+81/ 
data  ((b(i,J),i=1,13),J=7,9)/ 

&8 .  , 8 .  , 8 .  , 
& A •  , 8 .  ,8.  , 

&8.  ,8. 1444444444444446881, 8. 744444444444444e+88, 
&8. 31 6666666666666e+88, 8. 528243982439824e+88,8. 146341 46341 4634e881, 
&8.534878848788488e+88,8.  ,8.  , 

’!*  'I*  » 

&Q'.3  e+81  | Q*283333333333333e+Q1  1 8 1 548788487 88 4878e+88  \ 
&8. 731 7873178731 71 e88 1,8. 621 951 21 951 21 95e+88,Q.  , 
&8«  >8.  ,  8  «  , 
&8»  ,8.  ,8.  , 
&8.  ,8.  ,8.833333333333333e*81, 
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&8. 274398243982439e+89, 8. 731 7873178731 78e«81 ,8. 29121 951 21 951 22e+88/ 
data  ((b(i,j),is1f13),js1i,12)/ 

&8 .  >8.  >8.  , 

&8.  , 8 •  >8.  , 

&al  ,,el439824398243982e+88 ',9.1 46341 463414634e+99| 

&8 • 2926 82926 82926 8 e+89 ,9.  ,8.  , 

&8.  >8.  >8.  , 

&8 .  f 8.  >8.  , 

&8  a  ,  8  •  ,9a  , 

&8a  ,9.  ,9.  , 

&8  a  ,8a  ,9a  , 

&8  a  ,8a  ,8.  , 

&8 a  ,9a  ,9a  , 


6 8.  ,9a  ,1.8  e+88/ 

data  (ch(i) ,i=1 ,13)/ 

&  8  a  ,9a  ,8a  , 

69.  ,8.  ,9,32388 952389952 4 e+89 , 

&8.2571 428571 42857e+8i, 8. 2571 428571 42857e+88, 8. 321 428571 422857ef 81, 
&8.321 428571 428571 e«81 ,8.  ,8.488895238095238e*01 , 

&8.488995238995238et91/ 

data  erps , tol ,  to.lt/ 1 .0e*86 ,1.8ei03,1.0ef94/ 
t=t1 


18  do  28  i=  1 ,  n 
28  xdum(i)=x(i) 


y..  if (  iprint.eq.fl)  call  printr(x,t,tc) 

\j>  call  deriv(x,t,tc,f1) 

if ( iprint.eq.8 )  call  hamilt(t,x,f1 ,h) 
if(t.eqat2.or.dt.eq.8.8)  go  to  258 
38  continue 

b21=b(2, 1 )*dt 
do  42  i= 1 , n 

42  x(l)sb21»f1(i)+xdum(i) 
call  deriv(x,t , tc ,f2) 
b31=b(3,1 )#dt 
b32=b(3,2)«dt 

do  431=1, n 

43  x(i)=b31#f1(i)+b32*f2(i)+xdum( i) 
call  deriv(x,t,tc,f3) 
b41=b(4, 1 )#dt 
b43=b(4,3)*dt 

do  44  i=1 ,n 

44  x(i)=b41»f1(i)+b43*f3(i)+xdum(i) 
call  deriv(x,t,tc,f4) 
b51=b(5, 1 )*dt 
b53=b(5,3)*dt 

do  45  i= 1 , n 

45  x(i)=b51*f1(i)+b53#(f3(i)*f4(i))+xdum(i) 
call  deriv(x,t, tc ,f5) 

b61=b(6, 1 )*dt 
b64=b(6,4)*dt 
b65=b(6,5)*dt 
do  46  i=1 ,n 

46  X(i)rb61«f1(i)+b64»f4(i)+b65»f5(i)+xdum(i) 
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call  deriv(x,t,tc,f6) 
t>71  =b(7 , 1 ) 
b74=b(7,4)*dt 
b75=b(7,5)*dt 
b76=b(7,6)*dt 
do  47  i= 1  ,n 

47  X(i)=b71*f1(i)+b74*f4(i)+b75*f5(i)+b76*f6(i)+xdum(i) 
call  deriv(x,t,tc,f7) 

b81=b(8,1)*dt 
b85=b(8,5)*dt 
b86=b(8,6)»dt 
b87=b(8,7)*dt 
do  48  i=1,n 

48  X(i)=b87*f7(i)+b81*f1(i)+b86*f6(i)+b85*f5(i)+xdum(i) 
call  deriv( x ,t , tc ,f8) 

b91=b(9,1)*dt 
b94sb(9,4)*dt 
b95=b(9,5)*dt 
b96=b(9,6)*dt 
b97=b(9,7)*dt 
b98=b(9,8)»dt 
do  49  i= 1 , n 

49  X(i)=b97*f7(i)+b91*f1(i)+b98«f8(i)+b94*f4(i)+b96»f6(i)+ 
&b95*f5 ( i )+xdum( i ) 

call  deriv( x ,t ,tc ,f9 ) 
bl91sb(19, 1 )*dt 
b194=b(19,4)*dt 
b195=b(19,5)»dt 
b196=b(19,6)«dt 
b197=b(19,7)*dt 
b198=b(19,8)#dt 
do  5 9  i=1, n 

59  X(i)=«b32*f9(i)+b194»f4(i)+b197»f7(i)+b191«f1 (i)+b198*f8(i)+ 
4b196#f6 ( i )+b1 95*f5 ( i )+xdum( i ) 
call  deriv( x ,t , tc , fl 9 ) 

bl 1 1=b( 1 1 , 1 )*dt 

b114=b(11 ,4)*dt 
b115=b(11 ,5)*dt 
b116=b(1 1 ,6)*dt 
b117=b(11,7)*dt 
bl 1 8=b( 11,8 )*dt 
bl 1 9=bC1 1,9) «dt 
bill 9=b( 1 1 , 19 )*dt 
do  51  i= 1 , n 

f4ei)=b114»f4ei)+b115*f5(i) 

51  X(i)sb119*f9(i)+b1119*f19(i)+b117*f7(i)+b118*f8(i)+ 

&b1 1 1*f  1  ( i)+f4e i)+b1 16*f6( i)+xdum( i) 

call  deriv(x , t , tc  ,f 1 1 ) 
b121 sb( 12, 1 )#dt 
b126=b(12,6)»dt 

b128sbei2,8)*dt 
do  52  i=1 , n 

52  X(i)  =  b121*(f1(i)«f7ei))^b128»ef8ei)«r9(i))-t* 
&b126«(f6ei)9f19(i))+xdumei) 
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call  deriv(x,t,tc,f12) 
bl 31-b< 1 3, 1 )*dt 
b136=b(13,6)«dt 
b1j7=b(13,7)»dt 
b138=b(13,8)»dt 
b139=b(13,9)*dt 
b131«=b(13,1i)*dt 
do  53  1=1, n 

53  X(i)=b139*f9(i)+b1319*f18(i)+b131*f1(i)+b137*f7(i)+ 
&b138*f8( 1 )+f 12( i )*dt+f4( i )+b1 36#f6 ( i )+xdum( i) 
call  deriv(x,t,tc,f13) 
c6=ch(6)*dt 
c7=ch(7)*dt 
c9=ch(9)»dfc 
c12=ch(12)*dt 
do  199  1=1, n 

199  x(i)  =  xdum(i)  +  (e6#f6(i)+c7<t(f7(i)  +  f8(i))+c9,l(f9(i)  +  f'1S(i))  + 
&c12*(f12(i)+fl3(i))) 
dtlsdt 

if(t.lt.tc)  go  to  11« 

call  contrl (x , t , tc ,mu2,mu3 ,phi, cl ,pi) 

if(pi.ge,1,9.or.pi«le.99.5)  jflag=1 

1 f ( Jflag.eq*1.and.pi.ge.1.90or.pi.le.i8.5)  dt=dt1»9*98 

if (dt.eq. (9,98*dt1 ) )  go  to  39 

if ( jflag.eq.l .and, pi, le. 1 •fl.and.pl „ge .*9,5)  t2=t+dtl 
lie  go  to  299 

er=9*9 

do  159  J=1,n 

if(aba(x( j) ).lt.erps)  go  to  159 
es=abs((f1(j)+f11(j)9f12(j)if13(J))/x(j)) 
if (es.gt.er )erses 
159  continue 

er=er*abs(c12)+1 , 9o#29 
dt=tolt/er 
dt=dt1*dt*«9.125 
If ( er.gt,tol)  go  to  39 
if(abs(dt).gt.9,25)  dt=99,25 
299  t=t+dt1 

if(abs(dt),gt.abs(t2«t))  dt=t2«t 
go  to  19 
259  return 
end 


uj 
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subroutine  printr(x,t,tc) 
dimension  x(1§) 


call  contrl(x,t,te,niu2,mu3fphi,cl,pi) 
fpa=x(6)*57. 29577951 
psi-x(7)*57. 29577951 
b an k=ph 1*57,29577951 


write(5,1fl)  t , ( x( i) , i=1 ,5) f fpa,psi , (x( i) , is8, 1i) , bank, cl , pi 
19  format ( lx ,f7. 4,/, 2(5 (1x,e14.7> ,/),3(1x,f 14.7 ) ,2(1x,e14.7) ) 
return 


end 
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subroutine  hamilt(t,x,f ,h) 

dimension  x(10),f(10) 

hax(5)*f (1 )+x(8)*f (4)+x(2)*f (6)+x(3)*f (7) 

write(5,1>  x(5),f(1),x(8),f(4),x(2) ,f (6) ,x(3> ,f (7) ,h 

1  forma t(  lx,  ’hamiltonian '  ,3x,4(1xte14.7),/,5(1x,e14.7)) 
Esx(1)+x(4)»*2/2. 0/32.174 

write(5,2)  E 

2  format( 1 x , ’ speoific  energy  a  ',e14,7) 
return 

end 
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subroutine  deriv(x,t ,to,f) 
dimension  x(1i),f(1t) 
real  k,l,lh,mu2,mu3,nmax 

common  a,beta,clmax,cdo,g,k,nmax,tmax,w,xi4 

common  clnmax ,hf ifi , jflag ,qa , rosig 
call  contrl(x,t,tc,mu2,mu3»phi,cl,pi) 
lsqa#cl 
lhrl*beta 

daqa*(cdo+k*cl#*2) 
dh=d*beta 
co3g=cos(x(6 ) ) 
sing=sin( x ( 6 ) ) 
cossscos(x(7) ) 
sins=sin( x(7 ) ) 
f(1 )=x(4)*sing 

f (7)=g/x(4)/cosg*l/w*sin(phi) 

f(2)=cosg*(g*(x(8)+mu3*(1.Mbeta»x(4)*»2/2.8/g))|x(5)*x(4)) 

&«sing»x(2)»g/x(4)ix(3)*f(7)»tan(x(6)) 

f(3)=«.« 

f (4)=g*( ( tmax*pifd)/w*sing) 
f (5 )sg/w#(dh*( x(8 )+mu3)+lh*(hfifi/x(4)tmu2 ) ) 
f (6)=g/x(4)#(l/w#cos(phi)*cosg) 
f (8)=2.«*f(5)/x(4)/beta+(x(2)*f (6)+x(3)*f(7) )/x(4) 

&|sing*( x(5)+mu3*beta*x(4 ) ) 
f (9 )=x(4)#cosg*coss 
f (1fl)=x(4)*cosg*sins 
return 
end 
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subroutine  contrl(x ,t , to ,mu2,mu3,phi ,ol ,pi) 
dimension  x(1fi) 

real  k , 1 , nmax ,mu2,mu3  . „ 

common  afbeta^clmaXtCdo^gjkfnmax ,tmax,w,xi4 

common  clnmax,hfifi. Jflag,qa.rosig 

hfifiasqrt ( x(2)**2+( x(3)/oos( x(b) ) )**2 ) 

rosiga«.f923769*exp(beta*x(D) 

qa=8.5*rosig*x(4)**2*a 

phisatan2(«(x(3)/cos(x(6))) ,»x(2)) 

clnmaxs2.9*w*nmax/(rosig,,a*x(4)**2) 

if(x(8).lt.8.8)  cl=*hfifl/<2.8*k*x(4)»x(8)) 

if(x.(8)  .ge.8.8)  clsclmax 

if (cUgt.elmax)  clsclmax 

if (cl.gt.clnmax)  clsclnmax 

if(cl.lt.8.8)  clsfi.8 

dsqa*(cdo+k*cl*#2) 

mu2s8«8 

mu3=8«8 

if(t.lt.tc)  go  to  18  M  ,  v  v  , 

pis (d+w*sin(x(6) )*( 1 •88beta*x(4)**2/2*8/g) ) /tmax 
mu3=*x(8) 
go  to  38 

if (cl«eq#clnmax,and.t ,lt«tc.and.xi4.gt.x(4) ) 
&mu2=2.8#k*x(8)*clnmax+hfifi/x(4) 


go  to  28 
pisl.8 
pisM.5 


29 


if(x(8).eq.8.8) 
if(x(8).lt.8.«) 
if(x(8).gt.8.8) 
go  to  38 
if (xi4,gt,x(4) ) 
if (xi4,lt*x(4) ) 
&(1.8/x(4)+x(5)*sin(x(6))) 
if (f 18.gt.9.8)  pi= 1*8 
if(f1«.lt.8.«)  pi=«8.5 
38  return 
end 


f18=t (1.8/x(4)+x(5)*sin(x(6)) ) 
f18=2.8*g*hfifi*nmax/x(4)**2» 
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subroutine  ndmiss(x,xi,r ,rnorm) 
dimension  x(19) ,xi(1«) ,r(4) 
r(1)s(x(1)#xi(1)) 

r(2)s(x(6)#xi(6) ) *57*2957795 1*18*9 
r (3 )=(x(7)*xi (7)) *57.29577951*19.9 
r(4)s(x(4)txl(4) ) 
rnorms9.9 
do  19  i=1 ,4 

19  rnorm=rnorm+r(i)*r( i) 
r norms sqrt(rnorm) 
return 
end 


AFIT/GAE/AA/78D- S 


subroutine  deltas < x,dx,tfs,tes, If lag, ii,iii) 

dimension  x(1*),dx(4> 

dx(1)ax(1)»1,ie«Q5 

dx(2)sx(2)*1.9e«84 

dx(3)=tfs*1.9e|«3 

if ( lflag#eq*2)  dx(3)=x(3)*1.fle*94 

dx(4)=tcs*1 ,8e|83 

if(ii,eq,1.or«ii,eq.2)  x(ii)=x(ii)+iii*dx(ii) 

if ( ii*eq.3*and.lflag#ne.2)  tfsatfs+iii^dxCS) 

if(ii.eq.3.and.lflag.eq.2>  x(3)=x(3)+iii*dx(3) 

if(ii,eq.4)  tcs=tcs+iii*dx(4) 

return 

end 
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subroutine  mtxinv(a,b,c,n) 
dimension  a(n,n) ,b(n,n) ,c(n,n) 
do  5  i=1 ,n 
do  5  J= 1 ,n 
5  c(i, j)  =  a(i,  J) 
do  Ifl  i=1,n 
do  Ifl  J=1,n 

Ifl  b(i, j)s«#fl 

do  15  i=1  ,n 
15  b(i,i)=1.« 

do  35  i=1,n 
cstor1sc(i , i) 
do  2fl  j=1  ,n 
c(if j)sc(i, J)/cstor1 
29  b(i, j)=b(i, j)/cstor1 
J  =  9 

25  j=j+1 

if( J.eq.i.and.i.eq.n)  go  to  49 
if(j.eq.i)  j=j+1 
cstor2=c( j , i) 
do  39  ks1,n 

c( j ,k)=c( j,k)«cstor2*c(i,k) 

39  b(j,k)=b(J ,k)flcstor2*b( i ,k) 
if(j.ne#n)  go  to  25 
35  continue 
49  return 
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subroutine  search (xf ,xi,tf ,tc ,dz,r ,rnorm, lflag , ipr in t) 
dimension  xf  ( If )  ,xi(  16)  ,x( It)  ,r  (4)  ,rp(4)  ,rm(4)  ,dz(4) 
iters® 
choices® 

rnormms9999999999, 
steplsl .6 

step2sfl. 618633989 

iteratsf 

mflagsi 

1  itersiter+1 

gsrsstep1*step2#*(  it  erf  1#iterat) 
write(5,#)  iter,gsr 
do  2  js1,2 

2  x( j)sxf( j)+gsr*dz( j) 

if ( lflag. eq.2)  x(3)=xf(3)+gsr«dz(3) 
do  3  j=5,16 

3  x(j)sxf(j) 
tfsstf 

if ( lflag.eq, 1 .or ,lflag,eq,3 )  tfsstf+gsr#dz(3 ) 

tcss tc+gsr#dz(4 ) 

if ( tcs.gt.tfs)  tcsstfs 

if (tcs.lt. 6.6)  tcs=6.8 

call  assign(x ,xi> tf , lflag) 

call  integ(18,x,tfs,tcs,1 ) 

call  ndmiss(x ,xi,rp, rnormp) 

if (mflag.eq. 1 )  go  to  4 

if( rnormp.lt. ( 1 .25 *r norm) )  mflagsi 

if (mflag.eq, 6 )  go  to  4 

step1=gsr 

step2s®. 618633989 

iteratsiterfl 

4  write(5,5)  tfs , tcs , rp , rnormp 

5  format (lx, 'tfs, tcs, rp, rnormp' ,/,7(1x,f18.4)) 
if ( rnormp. gt.rnormm)  go  to  7 

do  6  is  1 , 4 

6  rm(i)=rp(i) 

if ( rnormp. lt.r norm)  rnormms rnormp 

7  tcsave=tcs 

if ( iter ,lt.26)  go  to  1 
tcs tcsave 

if(tc„gt„tf)  tcstf 
if (tc.lt, 6.661 )  tcs®. 811 
call  assign(xf ,xi ,tf , lflag) 
do  1/8  is  1,4 
176  r(i)srm(i) 
r norms r nor mm 
return 
end 
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The  repor£*deals  with  the  tradeoff  between  time-to-turn  and  the  change  in  specific 
energy  during  a  turn  of  180®,  for  an  aircraft  of  high  specific  power.  This  type  of 
aircraft  possesses  the  capability  to  sustain  flight  at  the  corner  velocity  where  both 
maximum  lift  coefficient  and  maximum  load  factor  occur  simultaneously.  However,  the 
classical  necessary  conditions  breakdown  on  this  corner  velocity  arc  and  an  additional 
constraint  must  be  defined  to  determine  the  optimal  control  histories.  The  report 
first  defines  the  necessary  conditions  for  a  generic  optimal  control  formulation  and 
then  applies  the  formulation  to  this  high  specific  power  aircraft  problem.  The  result 
is  a  three-point  boundary  value  problem  with  a  discontinuous  interior  corner  condition 
at  the  beginning  of  the  sustained  corner  velocity  arc.  All  state  and  costate  derivatives 
and  end  conditions  are  presented  with  numerical  methods  for  determining  both  minimum 
time  and  maximum  energy  gradient  solutions.  Anticipated  results  are  provided  for  the 
formulations  beginning  at  initial  conditions  both  above  and  below  the  corner  velocity. 
Results  are  presented  for  the  minimum  time  solution  only.  Recommends  ticnr.  for  further 


work  in  this  same  area  are  also  provided.  Q-'..'.  .r, 
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